Prior to explaining the invention, it is helpful to understand first the prior art of conventional Kalman recursions as well as the Fast Kalman Filtering (FKF) method for both calibrating a sensor system PCT/FI90/00122 (WO 90/13794) and controlling a large dynamic system PCT/FI93/00192 (WO 93/22625).
The underlying Markov (finite memory) process is described by the equations from (1) to (3). The first equation tells how a measurement vector y.sub.t depends on a state vector s.sub.t at time point t, (t=0,1,2 . . . ). This is the linearized Measurement (or observation) equation: EQU y.sub.t =H.sub.t s.sub.t +e.sub.t (1)
Matrix H.sub.t is the design (Jacobian) matrix that stems from the partial derivatives of actual physical dependencies. The second equation describes the time evolution of the overall system and it is known as the linearized System (or state) equation: EQU s.sub.t =A.sub.t s.sub.t-1 +B.sub.t u.sub.t-1 +a.sub.t (2)
Matrix A.sub.t is the state transition (Jacobian) matrix and B.sub.t is the control gain (Jacobian) matrix. Equation (2) tells how present state s.sub.t of the overall system develops from previous state s.sub.t-1, control/external forcing u.sub.t-1 and random error a.sub.t effects. When measurement errors e.sub.t and system errors a.sub.t are neither auto- (i.e. white noise) nor cross-correlated and are given by the following covariance matrices: EQU Q.sub.t =Cov(a.sub.t)=E(a.sub.t a.sub.t ')
and EQU R.sub.t =Cov(e.sub.t)=E(e.sub.t e.sub.t ') (3)
then the famous Kalman forward recursion formulae from (4) to (6) give us Best Linear Unbiased Estimate (BLUE) s.sub.t of present state s.sub.t as follows: EQU s.sub.t =A.sub.t s.sub.t-1 +B.sub.t u.sub.t-1 +K.sub.t {y.sub.t -H.sub.t (A.sub.t s.sub.t-1 +B.sub.t u.sub.t-1)} (4)
and the covariance matrix of its estimation errors as follows: EQU P.sub.t =Cov(s.sub.t)=E{(s.sub.t -s.sub.t)(s.sub.t -s.sub.t)'}=A.sub.t P.sub.t-1 A'.sub.t +Q.sub.t -K.sub.t H.sub.t (A.sub.t P.sub.t-1 A'.sub.t +Q.sub.t) (5)
where the Kalman gain matrix K.sub.t is defined by EQU K.sub.t =(A.sub.t P.sub.t-1 A'.sub.t +Q.sub.t)H'.sub.t {H.sub.t (A.sub.t P.sub.t-1 A'.sub.t +Q.sub.t)H'.sub.t +R.sub.t }.sup.-1 (6)
This recursive linear solution is (locally) optimal. The stability of the Kalman Filter (KF) requires that the observability and controlability conditions must also be satisfied (Kalman, 1960). However, equation (6) too often requires an overly large matrix to be inverted. Number n of the rows (and columns) of the matrix is as large as there are elements in measurement vector y.sub.t. A large n is needed for making the observability and controlability conditions satisfied. This is the problem sorted out by the discoveries reported here and in PCT/FI90/00122 and PCT/FI93/00192.
The following modified form of the State equation has been introduced EQU A.sub.t s.sub.t-1 +B.sub.t u.sub.t-1 =Is.sub.t +A.sub.t (s.sub.t-1 -s.sub.t-1)-a.sub.t (7)
and combined with the Measurement equation (1) in order to obtain the so-called Augmented Model: ##EQU1##
The state parameters can be estimated by using the well-known solution of a Regression Analysis problem as follows: EQU s.sub.t =(Z'.sub.t V.sub.t.sup.-1 Z.sub.t).sup.-1 Z'.sub.t V.sub.t.sup.-1 z.sub.t (9)
The result is algebraically equivalent to use of the Kalman recursions but not numerically (see e.g. Harvey, 1981: "Time Series Models", Philip Allan Publishers Ltd, Oxford, UK, pp. 101-119). The dimension of the matrix to be inverted in equation (9) is now the number (=m) of elements in state vector s.sub.t. Harvey's approach is fundamental to all different variations of the Fast Kalman Filtering (FKF) method.
An initialization or temporary training of any large Kalman Filter (KF), in order to make the observability condition satisfied, can be done by Lange's High-pass Filter (Lange, 1988). It exploits an analytical sparse-matrix inversion formula for solving regression models with the following so-called Canonical Block-angular matrix structure: ##EQU2##
This is the matrix representation of the Measurement equation of e.g. an entire windfinding intercomparison experiment. The vectors b.sub.1,b.sub.2, . . . ,b.sub.K typically refer to consecutive position coordinates e.g. of a weather balloon but may also contain those calibration parameters that have a significant time or space variation. The vector c refers to the other calibration parameters that are constant over the sampling period.
For all large multiple sensor systems their design matrices H.sub.t are sparse. Thus, one can do in one way or another the same sort of ##EQU3##
where
c.sub.t typically represents calibration parameters at time t PA1 b.sub.t,k all other state parameters in the time and/or space volume PA1 A state transition matrix (block-diagonal) at time t PA1 B matrix (block-diagonal) for state-independent effects u.sub.t at time t. PA1 c.sub.t typically represents "calibration" parameters at time t PA1 b.sub.t,k values of atmospheric parameters at gridpoint k (k=1, . . . K) PA1 A state transition matrix at time t (submatrices A.sub.1, . . . ,A.sub.K,A.sub.c) PA1 B control gain matrix (submatrices B.sub.1, . . . ,B.sub.K,B.sub.c).
If the partitioning is not obvious one may try to do it automatically by using a specific algorithm that converts every sparse linear system into the Canonical Block-angular form (Weil and Kettler, 1971: "Rearranging Matrices to Block-angular Form for Decomposition (and other) Algorithms", Management Science, Vol. 18, No. 1, Semptember 1971, pages 98-107). The covariance matrix of random errors e.sub.t may, however, loose something of its original and simple diagonality.
Consequently, gigantic Regression Analysis problems were faced as follows:
Augmented Model for a space volume case: e.g. for the data of a balloon tracking experiment with K consecutive balloon positions: ##EQU4##
Augmented Model for a moving time volume: (e.g. for "whitening" an observed "innovations" sequence of residuals e.sub.t for a moving sample of length L): ##EQU5##
Please note that the latter matrix equation has a "nested" Block-Angular structure. There are two types of "calibration" parameters c.sub.t and C.sub.t. The first set of these parameters, c.sub.t, can vary from one time to another. The second type, C.sub.t, of these parameters have constant values (approximately at least) over long moving time windows of length L. The latter ones, C.sub.t, make the Kalman filtering process an adaptive one. The solving of the latter parameters with the conventional Kalman recursions from (4) to (6) causes an observability problem as for computational reasons length L must be short. But with the FKF formulae of PCT/FI90/00122, the sample size can be so large that no initialization (or training) may be needed at all.
Prior to explaining the method of PCT/FI93/00192, it will be helpful to first understand some prior art of the Kalman Filtering (KF) theory exploited in experimental Numerical Weather Prediction (NWP) systems. As previously, they also make use of equation (1): EQU Measurement Equation: y.sub.t =H.sub.t s.sub.t +e.sub.t (linearized regression)
where state vector s.sub.t describes the state of the atmosphere at time t. Now, s.sub.t usually represents all gridpoint values of atmospheric variables e.g. the geopotential heights (actually, their departure values from the actual values estimated by some means) of different pressure levels.
The dynamics of the atmosphere is governed by a well-known set of partial differential equations ("primitive" equations). By using e.g. the tangent linear approximation of the NWP model the following expression of equation (2) is obtained for the time evolution of state parameters s.sub.t (actually, their departure values from a "trajectory" in the space of state parameters generated with the nonlinear NWP model) at a time step: EQU State Equation: s.sub.t =As.sub.t-1 +Bu.sub.t-1 +a.sub.t (the discretized dyn-stoch.model)
The four-dimensional data assimilation results (s.sub.t) and the NWP forecasts (s.sub.t) respectively, are obtained from the Kalman Filter system as follows: EQU s.sub.t =s.sub.t +K.sub.t (y.sub.t -H.sub.t s.sub.t) EQU s.sub.t =As.sub.t-1 +Bu.sub.t-1 (12)
where EQU P.sub.t =Cov(s.sub.t)=A Cov(s.sub.t-1)A'+Q.sub.t (prediction accuracy) EQU Q.sub.t =Cov(a.sub.t)=Ea.sub.t a'.sub.t (system noise) EQU R.sub.t =Cov(e.sub.t)=Ee.sub.t e'.sub.t (measurement noise)
and the crucial Updating computations are performed with the following Kalman Recursion: EQU K.sub.t =P.sub.t H'.sub.t (H.sub.t P.sub.t H'.sub.t +R.sub.t).sup.-1 (Kalman Gain matrix) EQU Cov(s.sub.t)=P.sub.t -K.sub.t H.sub.t P.sub.t (estimation accuracy).
The matrix inversion needed here for the computation of the Kalman Gain matrix is overwhelmingly difficult to compute for a real NWP system because the data assimilation system must be be able to digest several million data elements at a time. Dr. T. Gal-Chen reported on this problem in 1988 as follows: "There is hope that the developments of massively parallel supercomputers (e.g., 1000 desktop CRAYs working in tandem) could result in algorithms much closer to optimal . . . ", see "Report of the Critical Review Panel--Lower Tropospheric Profiling Symposium: Needs and Technologies", Bulletin of the American Meteorological Society, Vol. 71, No. 5, May 1990, page 684.
The method of PCT/FI93/00192 exploits the Augmented Model approach from equation (8): ##EQU6##
The following two sets of equations are obtained for Updating purposes: ##EQU7##
where, EQU s.sub.t =As.sub.t-1 +Bu.sub.t-1 (NWP "forecasting") EQU P.sub.t =Cov(s.sub.t)=A Cov(s.sub.t-1)A'+Q.sub.t (15)
but instead of EQU K.sub.t =P.sub.t H'.sub.t (H.sub.t P.sub.t H'.sub.t +R.sub.t) (Kalman Gain matrix)
the FKF method of PCT/FI93/00192 takes EQU K.sub.t =Cov(s.sub.t)H'.sub.t R.sub.t.sup.-1 (16)
The Augmented Model approach is superior to the use of the conventional Kalman recursions for a large vector of input data y.sub.t because the computation of the Kalman Gain matrix K.sub.t requires the huge matrix inversion when Cov(s.sub.t) is unknown. Both methods are algebraically and statistically equivalent but certainly not numerically.
However, the Augmented Model formulae are still too difficult to be handled numerically. This is, firstly, because state vector s.sub.t consists a large amount (=m) of gridpoint data for a realistic representation of the atmosphere. Secondly, there are many other state parameters that must be included in the state vector for a realistic NWP system. These are first of all related to systematic (calibration) errors of observing systems as well as to the so-called physical parameterization schemes of small scale atmospheric processes.
The calibration problems are overcome in PCT/FI93/00192 by using the method of decoupling states. It is done by performing the following Partitioning: ##EQU8##
where
Consequently, the following gigantic Regression Analysis problem was faced: ##EQU9##
The Fast Kalman Filter (FKF) formulae for the recursion step at any time point t were as follows: ##EQU10##
where, for k=1,2, . . . ,K, EQU R.sub.t,k =V.sub.t,k.sup.-1 {I-X.sub.t,k {X'.sub.t,k V.sub.t,k.sup.-1 X.sub.t,k }.sup.-1 X'.sub.t,k V.sub.t,k.sup.-1 }
 ##EQU11##
and, i.e. for k=0, EQU R.sub.t,0 =V.sub.t,0.sup.-1 EQU V.sub.t,0 =Cov{A.sub.c (s.sub.t-1 -s.sub.t-1)-a.sub.c.sub..sub.t } EQU Y.sub.t,0 =A.sub.c s.sub.t-1 +B.sub.c u.sub.t-1 EQU G.sub.t,0 =I.
The data assimilation accuracies were obtained from equation (20) as follows: ##EQU12##
Kalman Filter (KF) studies have also been reported e.g. by Stephen E. Cohn and David F. Parrish (1991): "The Behavior of Forecast Error Covariances for a Kalman Filter in Two Dimensions", Monthly Weather Review of the American Meteorological Society, Vol. 119, pp. 1757-1785. However, the ideal Kalman Filter systems described in all such reports is still out of reach for Four Dimensions (i.e. space and time). A reliable estimation and inversion of the error covariance matrix of the state parameters is required as Dr. Heikki Jarvinen of the European Centre for Medium-range Weather Forecasts (ECMWF) states: "In meteorology, the dimension (=m) of the state parameter vector s.sub.t may be 100,000-10,000,000. This makes it impossible in practice to exactly handle the error covariance matrix.", see "Meteorological Data Assimilation as a Variational Problem", Report No. 43 (1995), Department of Meteorology, University of Helsinki, page 10. Dr. Adrian Simmons of ECMWF confirmes that "the basic approach of Kalman Filtering is well established theoretically, but the computational requirement renders a full implementation intractable.", see ECMWF Newsletter Number 69 (Spring 1995), page 12.
The Fast Kalman Filtering (FKF) formulas known from PCT/FI90/00122 and PCT/FI93/00192 make use of the assumption that error covariance matrix V.sub.t in equations (9) and (13), respectively, is block-diagonal. Please see the FKF formula (19) where these diagonal blocks are written out as: ##EQU13##
It is clear especially for the case of adaptive Kalman Filtering (and the 4-dimensional data-assimilation) that the estimates of consecutive state parameter vectors s.sub.t-1, s.sub.t-2, s.sub.t-3, . . . are cross- and auto-correlated.
There exists a need for exploiting the principles of the Fast Kalman Filtering (FKF) method for adaptive Kalman Filtering (AKF) with equal or better computational speed, reliability, accuracy, and cost benefits than other Kalman Filtering methods can do. The invented method for an exact way of handling the error covariances will be disclosed herein.